﻿/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * Geophysical Computational Tools & Library (GCTL)
 *
 * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
 *
 * GCTL is distributed under a dual licensing scheme. You can redistribute 
 * it and/or modify it under the terms of the GNU Lesser General Public 
 * License as published by the Free Software Foundation, either version 2 
 * of the License, or (at your option) any later version. You should have 
 * received a copy of the GNU Lesser General Public License along with this 
 * program. If not, see <http://www.gnu.org/licenses/>.
 * 
 * If the terms and conditions of the LGPL v.2. would prevent you from using 
 * the GCTL, please consider the option to obtain a commercial license for a 
 * fee. These licenses are offered by the GCTL's original author. As a rule, 
 * licenses are provided "as-is", unlimited in time for a one time fee. Please 
 * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
 * to include some description of your company and the realm of its activities. 
 * Also add information on how to contact you by electronic and paper mail.
 ******************************************************/

#ifndef _GCTL_GEOMETRY3D_H
#define _GCTL_GEOMETRY3D_H

// library's head file
#include "../core/array.h"
#include "../core/matrix.h"
#include "../algorithm.h"

#include "tensor.h"
#include "node.h"
#include "edge.h"
#include "tetrahedron.h"
#include "triangle.h"
#include "block.h"
#include "tri_cone.h"
#include "sphere.h"
#include "tesseroid.h"
#include "prism.h"

namespace gctl
{
	namespace geometry3d
	{
		/**
		 * @brief      两个矢量以坐标原点为顶点的夹角（弧度）
		 *
		 * @param[in]  a     三维空间内的一个向量或实点a。
		 * @param[in]  b     三维空间内的一个向量或实点b。
		 *
		 * @return     两个矢量之间的夹角。
		 */
		double angle(const point3dc &a, const point3dc &b);

		/**
		 * @brief 计算三角形的面积
		 * 
		 * @param a 三维空间内的一个实点a
		 * @param b 三维空间内的一个实点b
		 * @param c 三维空间内的一个实点c
		 * @return double 面积
		 */
		double triangle_area(const point3dc &a, const point3dc &b, const point3dc &c);

		/**
		 * @brief      点到线段的距离
		 *
		 * @param[in]  line_start  线段的起点
		 * @param[in]  line_end    线段的终点
		 * @param[in]  dot         线段外的一个点
		 *
		 * @return     点到线段的距离
		 */
		double dot2line(const point3dc &line_start, const point3dc &line_end, const point3dc &dot);
		/**
		 * @brief      点在直线上的投影
		 *
		 * @param[in]  line_start  线段的起点
		 * @param[in]  line_end    线段的终点
		 * @param[in]  dot         线段外的一个点
		 *
		 * @return     点在直线上的投影
		 */
		point3dc dot_on_line(const point3dc &line_start, const point3dc &line_end, const point3dc &dot);
		/**
		 * @brief       点到平面的距离。
		 *
		 * @param[in]  c     平面上的一个点
		 * @param[in]  n     平面的单位外法线矢量
		 * @param[in]  d     平面外的一点
		 *
		 * @return     平面外的点到平面的距离
		 */
		double dot2plane(const point3dc &c, const point3dc &n, const point3dc &d);
		/**
		 * @brief Calculate the coefficients of a 3D plane using coordinates of three points on the plane. 
		 * The plane's equation is given as Ax + By + Cz + D = 0
		 * 
		 * @param x1  x coordinate of the first point
		 * @param x2  x coordinate of the second point
		 * @param x3  x coordinate of the third point
		 * @param y1  y coordinate of the first point
		 * @param y2  y coordinate of the second point
		 * @param y3  y coordinate of the third point
		 * @param z1  z coordinate of the first point
		 * @param z2  z coordinate of the second point
		 * @param z3  z coordinate of the third point
		 * @param A   coefficient A
		 * @param B   coefficient B
		 * @param C   coefficient C
		 * @param D   coefficient D
		 */
		void get_plane_coeff(double x1, double x2, double x3, double y1, double y2, double y3, 
			double z1, double z2, double z3, double &A, double &B, double &C, double &D);

		/**
		 * @brief Calculate the coefficients of a 3D plane defined by two points on a sphercial surface and the origin. 
		 * The plane's equation is given as Ax + By + Cz + D = 0
		 * 
		 * @param p1 point 1 on the spherical surface.
		 * @param p2 point 2 on the spherical surface.
		 * @param A  coefficient A
		 * @param B  coefficient B
		 * @param C  coefficient C
		 * @param D  coefficient D
		 */
		void get_plane_coeff(const point3ds &p1, const point3ds &p2, double &A, double &B, double &C, double &D);
		
		/**
		 * @brief      计算直线与平面的交点。
		 *
		 * @param[in]  face_p  平面上的一个点。
		 * @param[in]  face_nor  平面的单位外法线矢量。
		 * @param[in]  line_p  直线上的一个点。
		 * @param[in]  line_nor  直线的单位方向矢量。
		 *
		 * @return     返回的矢量。
		 */
		point3dc line_on_plane(const point3dc &face_p, const point3dc &face_nor, 
			const point3dc &line_p, const point3dc &line_nor);
		/**
		 * @brief      计算两个矢量之间的椭圆弧上的任意一点
		 * 
		 * 三维空间中任意两个以原点为起点矢量可视为一个三维椭圆弧上的两个点，给定这两个点与需要计算的第三点矢量
		 * 与起点矢量的夹角可以获得计算点坐标。
		 *
		 * @param[in]  v1    起点矢量
		 * @param[in]  v2    终点矢量
		 * @param[in]  phi   计算点矢量与起点矢量之间的夹角(弧度)
		 *
		 * @return     返回计算点坐标
		 */
		point3dc track_ellipse(const point3dc &v1, const point3dc &v2, double phi);
		/**
		 * @brief      点在平面上的投影
		 *
		 * @param[in]  face_dot  平面上一点
		 * @param[in]  face_nor  平面的法向量
		 * @param[in]  dot       平面外一点
		 *
		 * @return     平面外一点在平面上的投影点
		 */
		point3dc dot_on_plane(const point3dc &face_dot, const point3dc &face_nor, const point3dc &dot);
		/**
		 * @brief      距离平方反比形式的三角形内插值函数
		 *
		 * @param[in]  p     待插值点坐标
		 * @param[in]  p1    三角形顶点坐标1
		 * @param[in]  p2    三角形顶点坐标2
		 * @param[in]  p3    三角形顶点坐标3
		 * @param[in]  d1    顶点1的值
		 * @param[in]  d2    顶点2的值
		 * @param[in]  d3    顶点3的值
		 *
		 * @return     插值点的值
		 */
		double tri_interp_dis(const point3dc &p, const point3dc &p1, const point3dc &p2, 
			const point3dc &p3, double d1, double d2, double d3);
		/**
		 * @brief      球心角平方反比形式的三角形内插值函数
		 *
		 * @param[in]  p     待插值点坐标
		 * @param[in]  p1    三角形顶点坐标1
		 * @param[in]  p2    三角形顶点坐标2
		 * @param[in]  p3    三角形顶点坐标3
		 * @param[in]  d1    顶点1的值
		 * @param[in]  d2    顶点2的值
		 * @param[in]  d3    顶点3的值
		 *
		 * @return     插值点的值
		 */
		double tri_interp_ang(const point3dc &p, const point3dc &p1, const point3dc &p2, 
			const point3dc &p3, double d1, double d2, double d3);
		/**
		 * @brief      细分三角形面积比形式的三角形内插值函数
		 *
		 * @param[in]  p     待插值点坐标
		 * @param[in]  p1    三角形顶点坐标1
		 * @param[in]  p2    三角形顶点坐标2
		 * @param[in]  p3    三角形顶点坐标3
		 * @param[in]  d1    顶点1的值
		 * @param[in]  d2    顶点2的值
		 * @param[in]  d3    顶点3的值
		 *
		 * @return     插值点的值
		 */
		double tri_interp_area(const point3dc &p, const point3dc &p1, const point3dc &p2, 
			const point3dc &p3, double d1, double d2, double d3);
		/**
		 * @brief      距离平方反比形式的三角形内插值函数
		 *
		 * @param[in]  p     待插值点坐标
		 * @param[in]  p1    三角形顶点坐标1
		 * @param[in]  p2    三角形顶点坐标2
		 * @param[in]  p3    三角形顶点坐标3
		 * @param[in]  p4    三角形顶点坐标4
		 * @param[in]  d1    顶点1的值
		 * @param[in]  d2    顶点2的值
		 * @param[in]  d3    顶点3的值
		 * @param[in]  d4    顶点4的值
		 *
		 * @return     插值点的值
		 */
		double tet_interp_dis(const point3dc &p, const point3dc &p1, const point3dc &p2, 
			const point3dc &p3, const point3dc &p4, double d1, double d2, double d3, double d4);
		/**
		 * @brief      球坐标转直角坐标
		 *
		 * @param[in]  sp_ptr  球坐标点的指针
		 * @param      cp_ptr  直角坐标点的指针
		 *
		 * @return     执行是否成功
		 */
		int sph2car(const point3ds *sp_ptr, point3dc *cp_ptr);
		/**
		 * @brief      直角坐标转球坐标
		 *
		 * @param[in]  cp_ptr  直角坐标点的指针
		 * @param      sp_ptr  球坐标点的指针
		 *
		 * @return     执行是否成功
		 */
		int car2sph(const point3dc *cp_ptr, point3ds *sp_ptr);
		/**
		 * @brief     椭球坐标转直角坐标
		 * 
		 * @param lon    经度
		 * @param lat    维度
		 * @param r      短半径
		 * @param R      长半径
		 * @return point3dc  直角坐标系点   
		 */
		point3dc ellip2car(double lon, double lat, double r, double R);

		/**
		 * @brief      判断射线是否过三角形 在order_triangular_surface中调用
		 *
		 * @param[in]  origin        射线的起点
		 * @param[in]  direct        射线的方向
		 * @param[in]  tri_facet     待检测的三角形
		 * @param[in]  active_edge   长度为3的数组。表示对应三角形的边是否为激活的状态（激活表示如果射线与此边相交也算与三角形相交）
		 * @param[in]  cutoff_limit  计算叉乘的零值精度，小于此值即认为是0
		 * @param[in]  cross_loc     返回交点的位置
		 *
		 * @return     是否在三角形内
		 */
		bool crossed_tri_facet(const point3dc &origin, const point3dc &direct, const triangle &tri_facet, 
			const bool *active_edge, double cutoff_limit = 1e-10, point3dc *cross_loc = nullptr);

		/**
		 * @brief     监测点t1与t2是否在点p1，p2与p3构成的平面的同一侧
		 * 
		 * @param p1  平面上一点
		 * @param p2  平面上一点
		 * @param p3  平面上一点
		 * @param t1  待检测点
		 * @param t2  待检测点
		 * @return true  在同一侧
		 * @return false 不在同一侧
		 */
		bool same_side(const point3dc &p1, const point3dc &p2, const point3dc &p3, const point3dc &t1, const point3dc &t2);
		/**
		 * @brief      判断一个点是否在四面体内
		 * 
		 * 注意点落在四面体的顶点，边或者面上也算是在四面体内
		 *
		 * @param      test_p  待检测点
		 * @param      p1      第一个点
		 * @param      p2      第二个点
		 * @param      p3      第三个点
		 * @param      p4      第四个点
		 *
		 * @return     是否在四面体内
		 */
		bool node_in_tetrahedron(const point3dc &p1, const point3dc &p2, 
			const point3dc &p3, const point3dc &p4, const point3dc &test_p);
		/**
		 * @brief      对多面体表面的三角形的顶点进行排序，保证所有三角形顶点都是逆时针排序的
		 *
		 * @param      poly_tri_ptr   Pointer of the triangle array
		 *
		 * @return     status of the function
		 */
		void order_triangular_surface(array<triangle> &poly_tri, double cutoff_limit = 1e-10);
		/**
		 * @brief      获取四面体网络中的公共面及公共面两侧的四面体指针
		 *
		 * @param[in]  ele            四面体网络的元素数组
		 * @param      out_face       返回的公共三角形列表数组
		 * @param      out_neigh_ptr  返回的公共三角形两侧邻接四面体的指针数组的指针，
		 * 若指针不为空则返回邻接四面体指针数组。否则不计算此数组。
		 */
		void get_tetra_mesh_common_faces(const array<tetrahedron> &ele, array<triangle> &out_face, 
			array<tetrahedron*> *out_neigh_ptr = nullptr);
		/**
		 * @brief      计算四面体网络的顶点数量
		 * 
		 */
		size_t sort_node_number(const array<tetrahedron> &ele);
		/**
		 * @brief      计算四面体网络中与顶点相连的四面体的指针列表
		 *
		 * @param[in]  ele         四面体数组
		 * @param      neigh_list  返回的四面体指针列表
		 * @param      node_nump   四面体数组内顶点个数的指针
		 */
		void sort_node_neighbor(const array<tetrahedron> &ele, std::vector<std::vector<tetrahedron*> > &neigh_list, size_t *node_nump = nullptr);
		/**
		 * @brief      计算四面体网络内的各个四面体之间的相邻关系
		 *
		 * @param     ele    四面体数组
		 * @param      node_nump   四面体数组内顶点个数的指针
		 */
		void sort_tet_neighbor(array<tetrahedron> &ele, size_t *node_nump = nullptr);

		/**
		 * @brief      将随机位置的平面数据以距离平方反比的方式插值到输出节点位置上
		 * 
		 * @note       此函数的搜索区域形状为椭球 适合用于三维直角坐标系下的数据规则化或者切割剖面等
		 *
		 * @param      out_posi      输出顶点位置指针
		 * @param[in]  out_num       输出顶点个数
		 * @param      in_posi       输入数据的位置指针
		 * @param      in_val        输入数据的数值指针
		 * @param[in]  in_num        输入数据的个数
		 * @param[in]  search_xlen   搜索半径的x方向长度
		 * @param[in]  search_ylen   搜索半径的y方向长度
		 * @param[in]  search_zlen   搜索半径的z方向长度
		 *
		 * @return     规则节点位置的数值 排列顺序从左下上右上
		 */
		array<double> *p2p_dist(point3dc * out_posi, int out_num, point3dc *in_posi, double *in_val, 
			int in_num, double search_xlen, double search_ylen, double search_zlen);

		/**
		 * @brief      将随机位置的平面数据以距离平方反比的方式插值到输出节点位置上
		 * 
		 * @note       此函数的搜索区域形状为球面底的圆锥台 适合用于三维球坐标系下的数据规则化或者切割剖面等
		 *
		 * @param      out_posi      输出顶点位置指针
		 * @param[in]  out_num       输出顶点个数
		 * @param      in_posi       输入数据的位置指针
		 * @param      in_val        输入数据的数值指针
		 * @param[in]  in_num        输入数据的个数
		 * @param[in]  search_rlen   搜索区域的半径方向长度
		 * @param[in]  search_deg    搜索区域的球心角度大小
		 *
		 * @return     规则节点位置的数值 排列顺序从左下上右上
		 */
		array<double> *p2p_dist_sph(point3ds * out_posi, int out_num, point3ds *in_posi, double *in_val, 
			int in_num, double search_rlen, double search_deg);

		/**
		 * @brief      使用平面切割三角网络
		 * 
		 * @note       平面法线正方向的三角网络将被保留
		 *
		 * @param[in]  in_eles    输入的三维三角形数组（一般为成片的三角网）
		 * @param[in]  nor        平面的法线方向
		 * @param[in]  surf       平面上的一点
		 * @param      out_nodes  切割后的新网络顶点集
		 * @param      out_eles   切割后的新网络三角形元素集
		 */
		void cut_triangular_mesh(const array<triangle> &in_eles, const point3dc &nor, const point3dc &surf, 
			array<vertex3dc> &out_nodes, array<triangle> &out_eles);
	}
}
#endif // _GCTL_GEOMETRY3D_H